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Simulation of Non-stationary, Non-Gaussian Water 
Levels on the Great Lakes 

Todd L. Walton Jr. 1 ,M. ASCE and Leon E. Borgman 2 ,M. ASCE 

INTRODUCTION 

Non-stationary non-Gaussian time series having correlation structure are common 
to the coastal engineering field. Often one does not know much about the underlying 
process or processes causing the series to vary as it does and hence cannot predict 
its future values with any certainty. Future values of the series can run a very 
different course than past values and hence a direct application of past events to 
future scenerios is not always a safe course of action. To do proper engineering 
design within a probabilistic framework it is safer to extract as much pattern from 
the process as feasible including extremal behavior, and then to prepare various 
possible future scenerios of the process for design purposes by some type of stochastic 
simulation. It should be emphasized that these are not deterministic predictions of 
what will occur at a particular time in the future. Rather, the set of such scenerios 
represent typical future conditions which the design or operation must withstand. 
The design should be examined relative to each such scenario to establish a range 
of response from ”worst” to ’’best”. 

This paper will present one approach to simulation for a non-stationary non- 
Gaussian time series of hourly water levels on the Great Lakes. An earlier paper 

(Walton (1989)) discussed an autoregressive simulation approach in the time do- 

1 Hydraulic Engineer,U.S. Army Corps of Engineers, Coastal Engineering Research Center, 
U.S.A.E. Waterways Experiment Station, Vicksburg,MS. 

2 Professor, Statistics Department, University of Wyoming, Laraimc, Wyoming 
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main of the pseudocyclic behavior of the monthly mean water levels in the Great 
Lakes. One underlying need for dynamic simulation of water levels is to drive a cross 
shore sediment transport numerical model (see, for example, Kraus and Larson(1988) 
or Kriebel and Dean(1985)) for assessment of erosion volume versus frequency-of- 
occurrance curves. The dynamical aspects of the water level records are important 
to determination of such erosion volume/frequency-of-occurrance curves. As an ex¬ 
ample, a storm having a large peak water level but short duration can produce less 
erosion than a storm having a lower peak water level but longer duration. 

Examples of the processes to be simulated herein are provided in Figs.l and 2 which 
represent the average hourly water levels for the 1978 fall season (October through 
December) at two sites: Holland, Michigan; and Erie, Pennsylvania. These hourly 
water levels were measured at National Oceanic and Atmospheric Administration 
(NOAA) Great Lake water level gage sites, numbers 7031 and 3038, respectively, 
and were provided courtesy of the NOAA Great Lakes Acquisition Unit in Rockville, 
Maryland. 

It is obvious from the time series plots that the short term (hourly) lake level is 
not stationary as evidenced by the heteroscedastic nature of the series fluctuations 


about a changing mean level. This complex structure of water levels is due to 
the underlying non-stationary forcing functions of wind, barometric pressure, and 
precipitation, which drive the water levels from equilibrium. The complexity of the 
water level series is also due to the complex response of the Great Lakes basins to 
the forcing functions. One such basin response to extreme event forcing function is 


seiching wI 


cb is apparent in the water level records after an initial extreme event 


water level shock. 
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The non-Gaussian nature of these short term water levels is best seen in histograms 
of the water levels Figs. 3 and 4. These histograms depart from the Gaussian dis¬ 
tribution as might be expected considering the numerous complex forcing functions 
which drive the Great Lakes water levels. 

METHODOLOGY 

For simulation of the present water levels, it is apparent from the above dis¬ 
cussion that the simulation methodology utilized will have to deal with the non¬ 
stationary character of the data. In the present procedure this is accomplished via 
filtering/detrending/normal scores transformation procedures to reduce the data to 
stationarity, followed then by Gaussian correlated simulation, and, subsequently, by 
inverse operations to recapture the non-stationary non-Gaussian trending character¬ 
istics of the original data. The procedure deals with non-Gaussian data distribution 
functions of a varying unknown nature with tails representative of the measured 
extremes through a robust tail fitting procedure coupled with use of the empirical 
distribution function and a bootstrapping procedure for simulation purposes. 

As a basis for later comparison of the simulated data to the measured data, the 
following four plots are considered: the time series plot (shown in Figs. 1 and 2); 
the histogram ( Figs. 3 and 4); a spectral density plot, and an autocorrelation plot. 
Autocorrelation plots for the 1978 fall hourly water levels for the Holland and Erie 
sites are provided in Figs. 5 and 6 respectively. Spectral density plots for the same 
data as smoothed via a Gaussian window function are provided in Figs. 7 and 8. 

The entire procedure for simulating future scenerios of water levels consists of a 
pattern analysis procedure for reduction of data to stationary Gaussian correlated 
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random noise coupled with a simulation procedure for reconstructing future scenerios 
of water levels from simulated Gaussian random noise having the same underlying 
correlation structure as the original series. The two part analysis procedure is best 
described via block diagram as in Fig 9. A step by step discussion of the procedure 
follows. 


Stationarity Procedures 


The data is first low pass filtered to determine a time varying mean and time vary¬ 
ing standard deviation using a Gaussian shaped smoothing or convolving function 
of the form: 


w(t) = 


ex P( = if~) 


(i) 


where c is a constant to adjust the width of w(r). The Fourier transform of w(r) is 
given as: 


/ OO 

tu(r) exp (-i27r/r)5r 

■OO 


( 2 ) 


= exp(-7r/ 2 c 2 ) 


(3) 


hence tu(0) = 1/c and W(0) = 1. 

The effective width of w(t) is defined as the width of a square pulse which has 
height tu(0) and area equal to W(0) (i.e., the area under w(r)). Thus, the vaiue ”c” 
may be identified as the effective width of w(r). 

The filtered trend vtrend(t ) is defined in the continuous time domain as a convo¬ 
lution of the signal v(t) with the Gaussian smoothing function u>(r): 

f°° 

vtrend(t) = / w(r)v(t — r)br 
J— CO 


(4) 
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where w(r) is given a chosen effective width. Let Vtrend(f) be the Fourier transform 
of vt''end(t), and V(f) the Fourier transform of v(t). 

Vtrend(f) = F{vtrend(t)} = / v£rend(r)exp(—t'2w/r)5r (5) 

J — 00 

V(f) =F{v(t)} = f v(t) exp (—i2icfr)6r (6) 

*/ — CO 

Since convolution in the time domain is equivalent to multiplication in the fre¬ 
quency domain, we then have 

Vtrend(f) = W{f) • V(f) (7) 

It is convenient to make these operations in frequency domain to gain large savings 
in computer processing time possible through the use of the fast Fourier transform 
algorithm. 

In practice, only a finite piece of data is available and it >s natural to revert to 
discrete Fourier transform (DFT) definitions. In the discrete time domain the filtered 
trend vtrend(t) is defined as: 

N-i 

vtrend(n A t) — At ^ w(k A t)v((n - k) At) (8) 

k =o 

where the sequence v(n A (); n=0.N-l is made circularly periodic, with period N 

for times outside the defined data. That is; 

v((n — k) A t) = v((N + 7i — k) At) (9) 

u(NAi) = v(0) (10) 

This artificial periodicity effects ?i frrnrl(n A t) about, one effective width it the 
beginning and end of the tisne series. The one effective width arises because the 
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standard deviation of the Gaussian shape of w{r) is hence three standard 
deviations are approximately equal to one effective width, ”c”. The Gaussian curve 
approaches zero at three standard deviations i. jin the mean. 

Using discrete Fourier transforms with 

jv —1 

Vtrend(mA f) = At vtrend(n Af)exp( "* - ^ mn ) (11) 

n=o 

V{mAf) = At u(nA/)exp(-~—) (12) 

n=0 

W(m A/) = At w i n At)exp(-^p^) (13) 

n=0 

where vtrend{nAt), v{nAt) and tu(» At) are all made circularly periodic, the trend 
computation is made in the discrete frequency domain by the following equation: 

Vtrend(m A /) = W{m A /) • V(m A /) (14) 


where 

A 1 
} ~ (N A t) 

and the discrete time equivalent of Eq. 3 is 

W(m A f) = exp (-7r(m A /) 2 c 2 ) 


(15) 


(16) 


Letting the quantity EW be the number of A i increments in the effective width c, 
so that 

c = {EW) A t (17) 


it follows that 


W{m A f) == exp( 


—7r m 


?/ rprxrW 
(iS YY ) y 


N 2 


(18) 
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since 

(At)(A/) = l (19) 

Calculations are made very rapidly in the frequency domain using fast Fourier trans¬ 
form (FFT) techniques (Blahut,(1985)) and the resulting time series vtrend(n At); 
n=0,...,N-l is found by inverse fast Fourier transform (IFFT). In actual simulations 
the vtrend(n A t) series should be made longer than needed such that one effective 
width can be removed at the beginning and the end of the series due to the imposed 
circular periodicity assumptions as noted previously. 

The local variance can be defined in an analogous manner in the discrete time 
domain as a convolution of the instantaneous variance with the same smoothing 
filter w(t): 

N-i 

vvar(n A t) = w(k A t)ivar((n — k) At) (20) 

Jfc=0 

where the instantaneous variance is defined as 

ivar(n At) = (v(n At) — vtrend(n A t)) 2 (21) 

The local variance can then be found as before in the frequency domain by a mul¬ 
tiplication of the Fourier transform of w(n A t) by the Fourier transform of the 
instantaneous variance. A different effective width can be used for the variance 
filtering operation if so desired. On inverse Fourier transform, the time domain 
standard deviation series vsd(n A t); n=0,...,N-l can be found as 

vsd(n At) = \/vvar(n At) (22) 

Again as in the case of the '/trend series, in actual simulations the »>sd(n A i) series 
should be made longer than needed such that at least one effective width can be 
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discarded at the beginning and end of the series. In this regard the longer of the 

two effective widths (trend and variance filter effective width) should be used as the 

basis for discarding information in order to keep the two series properly aligned. 

The required stationary residuals are then found as; 

. a «(n A t) - vtrend(n A t) 

vresidfn At) = —- J -r. -- (23) 

v ' vsd(n At) v ' 

These residuals will locally have a mean of zero and standard deviation of unity. 

Normalization Procedure 

The residuals vresid(nAt) are then order ranked by size (smallest to largest) with 
the integer series rank(n A <);n=0,...,N-l being the rank number series where each 
value of rank(n A t) corresponds to the rank number of vresid(n A l). Hence the 
series rank(n A t);n=0,...,N-I ranges between 1 and N. An estimate of the distribu¬ 
tion function of the detrended and normalized (mean zero, variance one) data after 
ranking may be stated as: 

F(vresid(n A <)) = —— (24) 

Each vresid(n A t) is thus replaced by a standardized normal value ”zscore(n A f)” 
having the ranking associated with the empirical probability distribution fractile 
value This transformation is accomplished via the following equation: 

. zscore{n A t) = 1 ( ) (25) 

where ^ _1 () is the inverse normal distribution function which can be solved for 
numerically (Zelen and Severo(1964) eq. 26.2.23). 








Standardization Procedure 


As a final, usually very small adjustment, the zscore(n At); n=0,...,N-l series is 
then restandardized to a zero mean and unit variance via 

|([|A|)= ..m( , ,ai)-l.r (26) 

where zbar and zsd are the mean and standard deviation of the zscore series. 


Spectrum Calculation 


The ’’target” spectrum for simulation purposes is a smoothed version of the line 

spectrum of the series z(n A <);n=0.N-l. The discrete Fourier transform (DFT) 

of the z(n Af);n—0,...,N-1 series is: 

Z(m A /) = At y] z{n A t) exp ( *^ mn ) (27) 

n=0 


and the spectral line at / = is 


Sz(m A /) = 


Z{m A /) | 2 
N A t 


The spectral lines are smoothed over frequency in order to obtain adequate spec¬ 
tral density estimates. This is achieved with Gaussian smoothing in the frequency 
domain via the convolution operation: 

N-i exp(- rt fc ft/) a -) 

Sbar(m&f)=Y{ --} • Sz((m - k) A /) (29) 

* “ ^ C. / 


For numerical efficency this operation is actually performed in the lag time domain 


via the equation 


sbar(n A t) = sz(n A t) 


exp(-™^) 


(30) 
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where sz(n A t) is the inverse discrete Fourier transform (IDFT) of Sz(m A /), and 
EW is an effective width as before where EW = = cjN At. Returning to the 

frequency domain, the smoothed ’’target” spectrum to be used in the simulation 
process ,Sbar(m A /), is the DFT of sbar(n A t). 

Simulation Spectrum Generation 

The smoothed ’’target” spectrum is then standardized to unit area to require 

N -1 

£ Sbar(m A/) A/ = 1.0 (31) 

m=0 

Simulation is based on the random behavior of a real Gaussian covariance stationary 
periodic (period = N At) sequence zsim(n A <);n=0,...,N-l with a zero mean and 
variance equal to one. The DFT of such a random series is 

—i2irmn jy 

Zsim(m A f) = U m - iV m = At'^ / zsim(n A i)exp (———) for 0 < m < — (32) 

n=0 

where U m and V m are independent, normally distributed random variables with 
zero mean and variance given by Borgman(1973,1982,1990 in press) and Miller and 
Borgman(1985) as: 

T , . ... . ,, . ... . N A tSzsimfm A f) ... N .... 

Variance(U m ) = Varmnc?.{V m ) =- 2 - l / 0 < m < — (33) 

Variance(U m ) = N A tSzsim( 0) if m = 0 

N A f . N 
Variance(U m ) = N A tSzsim (—-—) if m = — 

N 

Variance[V m ) = 0 if m = 0, — 

where Szsim(m A f) is the desired line spectrum to simulate. For y < m < N, 


bJm — UN—m ond Vm — VN—m 


(34) 
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The frequency series Uo,Ui,Vi . Ux._ 1 ,Vu._ 1 ,Uu. is thus generated by generat¬ 

ing N independent standard Gaussian variates Z\ t Z 2t ....Zn and then calculating the 
frequency series in accordance with Szsim(m A f) = Sbar(m A /): 

U 0 = s/N A tSbar(0)Z 1 (35) 

t/i = s/N A tSbar{Af)/2Z 2 
V x = s/'N A tSbar(Af)/2Z 3 
U 2 = s/N A tSbar(2 A f)/2Z A 
V 2 = s/N A tSbar( 2 A f)/2Z 5 

etc. 

U*_i = ^NAtSbar((j-l)Af)/2Z N . 2 
VtL_ x - yN A tSbar((j - 1) A /)/2Z W -i 
Uz = s^N AtSbar{^ A f)Z N 

The simulated Gaussian time series with zero mean and unit variance zsim(n A 

t);n=0.N-l is then generated via the IDFT of the series Zsim(m A f) = U m - 

il/ m ;m=0,...,N-l with the equation: 

zsim(nAt) = Af ^(t/ m - ?V m )exp(^^) (36) 


Destandardization Procedure 

After generating the zero mean, unit variance time series zsim(n At);n=0,...,N-i, 
it is rescaled with the zscore(n A t);n=0,...,N-l series mean “zbar” and standard 
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deviation ”zsd”. This is done via the equation 

zscosim(n A t) = zsd • zsim(n At) + zbar (37) 

to create a simulated Gaussian time series zscosim(n A <);n=0.N-l having the 

same correlation structure, mean, and standard deviation, as the zscore(nA<);n=0,...,N' 
1 series. 


Distribution Function Transformation 

At this point the Gaussian simulated series zscosim(n A t);n=0.N-l is trans¬ 

formed back to the probability distribution function of the detrended series to obtain 
the resulting non-Gaussian trended simulated time series vsimres(n A<);n=0,...,N-l. 
This is done via back interpolation using the empirical distribution function in the 
central portion of the distribution along with a functional form for the tails of the 
distribution to make the distribution function more ’’robust”. The back interpola¬ 
tion from the simulated normal score time series zscosim(n A t) is achieved from 
inverting 

F(vsimres(n A t)) = <!>(zscosim(n A i)) (38) 

where $(z) is the distribution function for a standard normal probability law. This 
may be written as: 

vsimres(n At) = F~ 1 (^(zscosim(n A t))) (39) 

where jF _ 1 () is the inverse empirical probability distribution function. In practice, 
the step 


u = $(zscosim(n A i)) 


(40) 
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is performed with an algorithm from Zelen and Severo(1964) (eq. 26.2.17) and then 
the operation 

vsimres(n At) = F~ x (u) (41) 

is achieved by direct interpolation within the ranked detrended original data, at 
least for l/(iV + 1) < u < N/(N + 1). 

If u < 1/(JV+1) (lower tail) or u > N/(N + 1) (upper tail), another procedure was 
selected. It was found that the upper and lower tails of many common probability 
laws such as the normal, exponential, gamma, and lognormal could be fitted e'ther 
exactly or with high accuracy of approximation by a representation of the form: 

v = z\ + (a + bT) c for upper tail (42) 

or 

v = z\ - (a + bT) c for lower tail 


where 

T = \/-2/n(l - U) upper tail unbounded above (43) 

T = U - U\ upper tail bounded above 
T = \/-2 ln(U) lower tail unbounded below 
T — U\ —U lower tail bounded below 

The constants a,b,c, and z\ are based on extremal values of the F(vresid(n A t)) 
empirical distribution function estimate. 

For the upper tail unbounded above, let z\ < Z 2 < zz be selected values within the 
upper tail of the ranked vresid(nAt) values, and u\ <uz < uz be the corresponding 
values of rank/(N+l). 
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If the upper tail is bounded by some value z top , then set z 3 — z top and u 3 = 1.0. 
The other two sets of (z,u) are obtained as before. 

For the lower tail unbounded below, the procedure is very similar. Let z 3 < zi < zi 
be selected within the lower ranked values of vresid(n Ai) and W 3 < uj < ui be the 
rank/(N+l) values associated with them. 

If the lower tail is bounded below by zjoitom) then set z 3 =s and u 3 ~ 0.0. 

The other two pairs of (z, u) are picked as for the unbounded below case. 

The zi in the previous tail representation formula coincides exactly with the z\ 
values assigned for each tail here. The «i values needed for the bounded cases in 
each tail are the corresponding «i values selected here. 

This leaves the a,b,and c values which can be computed by: 


c 




(44) 


ln(b) = 


In 1 z 3 — zi 1 
c 


-In ( T 1 -T 1 | 


(45) 


o = - 6 T 1 


(46) 


where T\ is the value of T computed from uj in the previous definition of T. 

The direct empirical interpolation is performed for lowertailui <u< uppertailui, 
while the tail approximation formulas are used for values outside this interval. 

Although the upper and lower tail represent a very small fraction of the simulated 
values, they are important in many engineering problems where extremes are signif¬ 
icant. In most published ’’bootstrap” procedures, the simulations are restricted to 
lying between the maximum and minimum of the original data. In predicting the 
future, these methods appear artificial, since it seems obvious that extreme values 
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outside the historical data base may occur eventually unless there is some sort of 
finite bounds due to the physical constraints. 

The procedure presented represents a type of tail-extended bootstrap simulation 
generalized with a normal score transformation. A great deal of research and com¬ 
parison with common population probability laws was involved in selecting the tail 
formulation recommended above (Borgman,(1989)). It appears reasonable under 
circumstances examined to date. 

Nonstationarity Transformation 

The last transformation necessary to produce a simulated series vsim(nAt);n=0,...,N- 
1 is an inverse transformation using the computed vtrend(nAt) and vsd(nAt) series. 
The final simulated non-stationary, non-Gaussian time series is found as; 

vsim(n A i) = ~d(n A t) • vsimres(n At) + vtrend(n A t) (47) 

This operation is performed in the time domain. 

RESULTS 

Results of the simulation procedure for one simulation run for the Erie site is 
presented in time series plot, Fig. 10, histogram Fig. 11, autocorrelation Fig.12, and 
smoothed spectrum (via a Gaussian window) Fig. 13. These figures can be compared 
to the actual series plots in Figs.2,4,6, and 8. The low pass trend and standard 
deviation series for the Erie site are shown in Fig. 14 and Fig. 15 respectively. The 
effective filter width utilized for the trend series was 20 hours while the effective filter 
width utilized for the standard deviation series was 50 hours. As the effective width 
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chosen is somewhat of an ad hoc procedure, no justification for these choices of filter 
width is given. A natural extension of these procedures is to analyze and synthesize 
the smoothing trend and local standard deviation series in a similar manner rather 
than use the deterministic output directly from the data analysis. This would not 
involve any new theory but rather would introduce a hierarchy of simulations. It is 
not being done here because it would introduce unnecessary complexity in explaining 
the method. 

Results of the simulation procedure for three simulation runs at the Holland site 
are provided in Figs. 16 thru 19 for the time series, histogram, autocorrelation 
and smoothed spectra (smoothed via a Gaussian window). The low pass trend and 
standard deviation series for this site were computed with the same effective filter 
width as the Erie site. Results are shown in Figs. 20 and 21. 

CONCLUSIONS 

The results of the procedure seem to ’’mimic” the original data in a reasonable 
fashon while still incorporating process noise consistant with the original series. The 
method can analyze and simulate a large time series extremely fast. The method 
allows for exact duplication of the univariate probability law contained in the data 
behavior and is thus capable of producing non- gaussian process behavior related 
to skewness and higher order moments. It is understood that some of the higher 
order bivariate moments may not be properly reproduced by this transformation to 
gaussian (and the inverse transformation). Applications for which the higher order 
bivariate (trivariate, etc.) moments are important should consider other techniques, 
perhaps based on Volterra series. The method also allows for the treatment of 
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non-stationarity in mean and variance either through including the bursts of non- 
stationarity of data in deterministic fashion or through a hierarchy of simulations as 
suggested in the text. 
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APPENDIX II. NOT/TION 
a,b = empirically fit constants; 

c = Gaussian time domain weighting function constant, also empirically fit con¬ 
stant; 

c/= Gaussian frequency domain weighting function constant; 

EW = effective width of filter; 

/ = frequency; 

ivar(n A t)= discrete time instantaneous variance series; 
m, n = index counters; 

N = number of signal samples; 

rank(n At) = discrete time integer ranking series; 

sbar(n At) = inverse discrete Fourier transform of Sbar(m A f) series; 

Sbar(m A f) = smoothed spectral line function of z series; 

Sz(m A f) = spectral line function of z series; 

Szsim(m A /)=desired spectral line function series for simulation; 
t = time; 

T = tail function value; 

U{‘,i = l,N — fractile value 

U m = real part of Zsim(m A f) series; 

v(t) = signal function, time representation; 

V(f) = signal function, frequency representation; 

V m = imaginary part of Zsim(m A f) series; 

vresid(n At) = discrete time stationary signal series; 

vsd(n At) = discrete time standard deviation series; 

vsim(n At) = discrete time simulated signal series; 

vsimres(n At) = discrete time simulated stationary signal series; 

vtrend(t) = filtered trend function, time representation; 

vtrend(n At) = discrete time trend series; 

Wrend(/)=filtcred trend function,frequency representation; 

Vtrend(m A /)=discrete time Fourier transform of vtrend(n At) series; 
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vvar(n A t) = discrete time smoothed variance series; 
w(r) = weighting function, time domain representation; 
w(n A f)= discrete time weighting series; 

W(f) = weighting function, frequency domain representation; 

W{m A /)=discrete time Fourier transform of w(n A t) series; 
z(n A t) = discrete time standardized zscore series; 

Z(m A /)= discrete Fourier transform of z series; 

Zi\i = 1 ,N independent standard Gaussian random variable; 
zbar s= mean of the zscore series; 

zscore(n At) = discrete time standardized normal score series; 
zsd = standard deviation of the zscore series; 

zscosim(n At) = discrete time restandardized simulated Gaussian series; 
zsim(n A t)= discrete time simulated signal series; 

Zsim(m A f) = discrete time Fourier transform of zsim(n A t) series; 

F() = empirical distribution function; 

F~ 1 () = inverse empirical distribution function; 

T() = Fourier transform notation; 

$() = standard normal distribution function; 

$ _1 () = inverse normal distribution function; 

A/ = frequency step; 

At = time step; 
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Figure 1. Hourly Water Level, Oct. 1 - Dec. 31, 1978, Holland, HI. 
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Figure 2. Hourly Water Level, Oct. 1 - Dec. 31, 1978, Erie, PA. 
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Figure 3. Histogram of Hourly Water Levels, Oct. - Dec. 1978, Holland, MI. 
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Figure 4. Histogram of Hourly Water Levels, Oct. - Dec. 1978, Erie, 


PA. 
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Figure 8. Smoothed Spectra, Hourly Water Levels, Oct. - Dec. 1978, Erie, PA. 
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Figure 9. Block Diagram of Simulation Procedure. 
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Figure 10. Simulated Hourly Water Levels, Oct. - Dec., Erie 











Figure 11. Histogram of Simulated Hourly Water Levels, Oct. - Dec., 

Erie, PA. 
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Figure 12. Autocorrelation, Simulated Hourly Water Levels, 
Oct. - Dec., Erie, PA. 
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Figure 13. Smoothed Spectra, Simulated Hourly Water Levels, 
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Figure 14. Trend, Hourly Water Levels, Oct. - Dec. 1978, 
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Figure 15. Local Standard Deviation, Hourly Water Levels, 
Oct. - Dec. 1978, Erie, PA. 




580.98 

579.75 
579.50 

579.25 

579.99 

578.75 

578.59 

578.25 
578.09 

577.75 

577.59 


Holland 7931 fall '78 si*3 



1.00 


432.80 


864.60 1296.40 1728.20 2160.00 


HOUR (1 = 0100 OCT. 1, 1878) 


Figure 16. Simulated Hourly Water Levels, Oct. - Dec., (3 

Holland, MI. 


Runs), 




































Figure 18 
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Figure 19. Smoothed Spectra, Simulated Hourly Water Levels, 
Oct. - Dec., (3 Runs), Holland, MI. 
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Figure 20. Trend, Hourly Water Levels, Oct. - Dec. 1978, Holland, MI. 
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Figure 21. Local Standard Deviation, Hourly Water Levels, Oct. - Dec. 1978, 

Holland, MI. 















